NONLINEAR ANALYSIS OF THE SOLUTIONS OF THE 
HASEGAWA WAKATANI EQUATIONS 

LINDA STALS * 

Abstract. The Hasegawa— Wakatani models are used in the study of confinement of hot plasmas 
with externally imposed magnetic fields. 

The nonlinear terms in the Hasegawa- Wakatani models complicate the analysis of the system as 
they propagate local changes across the entire system. Centre manifold analysis allows us to project 
down onto much smaller systems that are more easily analysed. Qualitative information about the 
behaviour of the reduced system, such as whether it is stable or unstable, can be used to predict the 

^~^ behaviour of the original full system. We show how the simple structure of the linear part of the 

T~H Hasegawa-Wakatani equations can be used to define these projection operators. 

C 2 ^ The centre manifold analysis will be used on a few examples to highlight certain properties of 

fSj the Hasegawa-Wakatani models. 

'~{ Key words. Hasegawa-Wakatani equations, centre manifold theory, nonlinear ordinary differ- 

h"~5 ential equations. 

1. Introduction. A number of papers focus on interpreting the numerical simu- 
lations of the Hasegawa-Wakatani equations [TJ [31 [SI [12] , but few have tried to analyse 
the system's behaviour and those that have are applicable under limited conditions 
[H [ini [in [E] ■ In this paper we develop a centre manifold analysis that predicts and 

,^ explains many interesting aspects of the system's behaviour, including examples of 

^ unstable solutions. The particular structure of the Hasegawa-Wakatani system allows 

Cj us to explicitly write down an approximation to the centre manifold, without the need 

I of any expensive computation, and thus we can study the behaviour of the system 

over a wide range of parameters. See Section [4J 

The author has written a code to solve the Hasegawa-Wakatani equations and the 
initial motivation behind the work presented in this paper was the validation of that 
code. Sections [5] and [6] relate the numerical simulations of the Hasegawa-Wakatani 
T-H equations to the centre manifold analysis. 

2. The Hasegavifa— Wakatani model. The Hasegawa-Wakatani model 8^ was 
^^ designed to extend the one field Hasegawa-Mima model |7i. The equations in a two- 
^^ dimensional domain couple the flow field given by the electrostatic potential 4> with 
,__! the density p by 

.> |v20+[0,V2</.]=a(0-p)-(-imv2(P+i)0 (2.1) 






> 



i;H _p+[^,p]=«(0_p)_^_Z_(_l)P/?^V2V (2.2) 



^ -p+[<^,p]=a(0-p)-.- 

where = (j){x , y , t) , p — p{x,y,t), V^ — d^/dx^ + d^/dy"^ is the two-dimensional 
Laplacian and p = 1, 2 is the order of the dissipation operator. The Poisson bracket 
[., .] defined by 

\f 1 ^ ^^ _ ^^ 
dx dy dy dx ' 

gives the convective derivative. The parameters, a, jS^, jSp and k are usually taken to 
be non-negative, but to simplify the analysis we will assume a, (3^, /3p are all positive. 
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The coefRcients (3^ and /3p denote viscosity and diffusion coefficients respectively. 
These terms, often called hyper-diffusions, may be small but are necessary to damp 
fluctuation energy that reaches the smallest scales. Physically, the k term feeds energy 
to the system, while the energy is dissipated by the parallel resistivity (a) and hyper- 
diffusions. We refer the reader to [H [8l [TOl [12] for more discussions on the physical 
interpretation of the system. 

2.1. Spectral discretisation of the HasegawaWakatani system. It is ap- 



propriate to use a Fourier spectral method to discrctise Equations (2.1) and (2.2 1 
since much of the essential physics of interest is contained in the energy distributions 
over scales of the motions. 

We assume that (j) and p are periodic over the domain fJ — [— tt, vr] x [— tt, tt]. 
Following the structure given in [B] define 






where co^ = e'^, uiy = e'^, k — {kx,ky), fc^' = (fc^ + ky)'^ and k^, ky are integers. 
We assume that 4>,p G H^'^^. Given m = 2n for some positive integer n, let the 
corresponding discrete space H^:^^ C HP^^ be defined as 

HP+' in) ^L:ue HP+' {n),u= £ U^ u;^ w^-" i . (2.3) 

Consequently (j)„i G Hm and /9,„ G Hm may be written as 

(prn{x,y)= J2 ^fe^^Wy". (2-4) 

n-1 

Pm{x,y)= Yl RkUJ^ujy'- (2.5) 



Substituting expansions (2.4) and (2.5) into the Hasegawa-Wakatani equations 



(2.1) and (2.2), making use of the differentiation properties of the Fourier transforms 



and collecting like terms gives 

-fc2^$fe + (7V(*, DM^))k = « ($fc - i?fe) + fc2(''+i)/3^$fc, 

Q (2.6) 

-Rk + (A^(*, R))k - a ($fe - Rk) - ikyK $fc - k^PppRk, 

and Dat = diag(fc^). The nonlinear terms N{^, Dn^) and N{^, R) are given by the 
convolution sums; 

n-l 

{N{^,Dn^))^^ Y^ {k^qy-qxky){{k^-'qxf + {ky~qyf)^q^k-q, (2.7) 



and 



(iV(*, R))^ = Y^ {k^qy - q^ky)<PqRk- 



It is convenient to write down the nonlinear term as a convolution sum for the theo- 
retical analysis, but in practise these terms are evaluated using fast Fourier transfor- 
mations. 

To fit the standard definition of an ordinary differential equation we would like to 
divide both sides of Equation (2.6) by — fc^, but k^ is zero when fc = (0,0). In Fourier 



space the (0, 0) mode is equivalent to the (2n, 2n) mode, which suggest the following 
definition, 



ki = 



fe2 

8n2 



if fc V 
if fc2 = 



The Hasegawa-Wakatani Equations in Fourier space are 

ai*'= = fci fci ^^'^ '''=' 



(2.9) 



-Rk = -(iV(*, R))k + a ($fe - Rk) - ikyK^k - Ppk'PRk, (2.10) 



for — n < kx, ky < n. 



From the definition of i/,„, see (2.3 1, we require ^o — Rq — for all values of t. 



According to Equations (2.9) and (2.10), if the initial values of $o and Rq are they 
remain 0. 

2.2. The numerical solution of the HasegawaWakatani equations. The 



system of ordinary differential equations defined by (2.9) and (2.10) are stiff and 



appropriate numerical schemes must be used to extract long term behaviour patterns. 
The author has studied the applicability of several numerical schemes for the solution 
of the Hasegawa-Wakatani equations, see [T7], and has written a Fortran 90 code 
that solves the equations using a fourth order implicit-explicit variable time step 
BDF scheme [HI [HI [19] . The results reported in Section [s] were obtained through 
the use of this code. High frequency components are damped by using a 2/3 dealising 
scheme similar to the one described by Pedersen et. al. [T^]- Pedersen et. al. also 
use a semi-implicit scheme. 

Another code referenced by a number of papers [H [HH] has been developed by 
Scott [131 [H]- In UM Scott argues that implicit methods are not appropriate for 
the Hasegawa-Wakatani equations. Our analysis in [17] shows that the Hasegawa- 
Wakatani equations are stiff and the bound on the time step size required to ensure 
stability makes the use of explicit methods impractical. We have found the implicit- 
explicit BDF method to be reliable and efficient. 

3. Eigendecomposition. The idea of the centre manifold analysis is to project 



the large system defined by Equations (2.9) and (2.10) onto a much smaller system 



that retains much of the qualitative type of behaviour of the original system. In order 
to find such a projection we firstly need to look at the eigenvalues and eigenvectors 
of the linear part the Hasegawa-Wakatani equations. 



Rewrite Equations (2.9) and (2.10) as 



Ft 



R 



= L 



R 
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F{^,R) 



(3.1) 



where 



L = 



A B 
C D 



A^dmg[-^-l3^k^py B^diagl^ 



C = diag (a — ikyK) , D = diag (—a — (3pk'^^) 



and 



F{^,R) 



(diag(fe^))-i 0' 
/ 






Recall — n < k^, ky < n. So $ G 



Re 



and L E 



3.1. Eigenvalues of linear system. A linear approximation to ( 2.9 1 and ( 2.10 1 
about the point (0, 0) is given by 



Ft 



= L 



Given the block structure of the L matrix it is possible to explicitly write down 
its eigenvalues and eigenvectors. Such information is crucial to the centre manifold 
analysis as it shows which modes influence the long term behaviour of the system. 
Suppose that A is not an eigenvalue of A and consider 



A/I 



A -XI B 

C D-XI 

A- XI B 

{D-XI)-C{A-XI)-^B 



I e 



'™ is an identity matrix. The above determinant will be zero when 
{A - XI){D - XI) -CB = X?I - X{A + D) + AD-CB = 0. 
Hence the 2to^, not necessarily distinct, eigenvalues are given by 



A^ = 



{A + D)k ± V {{A + D)kr -4 {AD- CB)k 



(3.2) 



(3.3) 



The following lemma confirms that the matrix A — XI is non-singular. 

Lemma 3.1. Eigenvalues of L are not eigenvalues of A. 

Proof. Suppose one of the eigenvalues of L, say A, was an eigenvalue of A. Then, 



as A is a diagonal matrix, Af, — A = for some k. Substituting this into (3.2 1 gives 



CBh — 0, which intern implies a = 0. A contradiction on the assumption that a > 0. 
D 

For the stability analysis we are interested in the sign of the real part of the 
eigenvalues. Firstly note that 



iA + D)k = -, ,^2 



+ a + (/3^ + Pp)k^A e M < 0, 



so 3fi(A^ ) < for all k. For terms under the square root we see that 

= -^ \a[3pk^P + af3^k^^P+^^ + (3^Ppk^^^P+^^ + iakyK 
fe_l_ L 

Consequently, if k is large enough, 5R(A^) > for small values of k^ and 3fi(A^) < 
for large values of k^ . 

When k = (0,0), Aj = 0. We always assume that the initial values of *(o.o) and 
Rujfi) are and remain over time, avoiding the issue of non-uniqueness arising from 
the eigenvalue. 

For large values of fc^, 

,± -(/30 + /3p)fc^P±|/30-/3p|fc^P 
K ^ • 

In which case the eigenvalues are approximately —f3pk'^P or —(3^k'^P, suggesting that 
the high frequency components of the solution quickly approach zero. 

Studies of the linear part of the Hasegawa-Wakatani equation have already been 
carried out by several authors, see [3J HI [TOl [H], although in these references the 
analysis was not expressed explicitly in terms of eigenvalues. 

3.2. Eigenvectors. Once again, due to the block diagonal structure of the L 
matrix we can explicitly write down the eigenvectors. 

Consider one of the eigenvalues A^ (or A^) and let w e M'" be defined by 



Wi 



1, ifj = fe 

0, otherwise 



Now let u S C™ be given by 

u+ = -{A-X+)-'^Bw. 
(or 

u- = -iA-X^r'Bw.) 



By Lemma 3.1 (^ — A^ /) is nonsingular and as it is a diagonal matrix the inverse 
can be easily found. The eigenvector corresponding to the eigenvalue A^ is [m* w\ . 

4. Analysis of nonlinear behaviour. Section [3] focussed on the behaviour 
of the linear part of the Hasegawa-Wakatani equations, which is a necessary first 
step to understanding the nonlinear behaviour. The highly structured nature of the 
Hasegawa-Wakatani equations means we are able to systematically use centre mani- 
fold theory to predict the behaviour of the nonlinear system. 

To apply the ideas behind centre manifold theory we firstly need to build an 
appropriate suspended system. 
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4.1. Suspended system. The matrix in (3.1) depends on the parameters a, k 
and /3^ and /3p. Let e represent one of the previously mentioned parameters and set 
e = i — e* . We wiU use e* in Section [5] to reference a particular choice of parameters 
that gives eigenvalues with zero real parts. A suspended system corresponding to the 



3.11 is 












d 
Ft 


R 

di 


- =0. 


R 


+ eG 


* 
R 



F{^,R), 



(4.1) 



The G matrix represents the interaction of the parameter e with R and $. For 
example, if e = a, we can rewrite (|2.9[) and (|2.10[) as 



d_ _ (JV(#, JjV^))fc a*{^k-Rk) 
dt '' k\ k\ 

{a - a*) ($fc - Rk) 



Afc2P$, 



dt 



Rk = -(Af(*, R))k + a* ($fc - Rk) - ikyK<^k - Ppk'PRk 
+ {a - a*) {<^k - Rk) , 



and 



G 



-(diag(fc^))-i (diag(fe^))-i 
/ -/ 



The eigenvalues and eigenvectors of L{e*) are discussed in Sections 3.1 and 3.2 



To help with the notation we will place the eigenvectors into the columns of the matrix 
P. In particular, 



P 



-{A - diag(A+))-ii3 -{A - diag(A^ ))-iS 
/ / 



Due to the simple structure of P we can find 



-B-'{A - dmg{X^))E{A - diag(A+)) -B-\A - die.g{X+))EB 
B-\A - diag(A^ ))i?(A - diag(A+)) B-\A - dmg{X^))EB _ 



where 



i?=(diag(A+)-diag(A^)) \ 



which is non-singular. The entries of P ^ niay look complicated, but they are all 
diagonal matrices. 
Let 



A = 



diag(A^ 







diag(Afc) 

We then have the following eigendecomposition 

L{€*)^PAP-\ 
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Now let n be a permutation matrix such that the first a entries along the diagonal 
of A = nAn^^ have zero real parts and the remaining b — 2rn^ — a entries have non- 



zero real components. Multiplying both sides of (4.1) by Q = HP ^ gives 

d_ 
where 



-X = KxX + Fx{X,Y) 
-Y = KyY + Fy{X,Y), 



Q 



A 



r*i 




rxi 


R 




Y 


\^x 


■ 







A 


Y 


; 



XeCM'e 



Ax e 



,Ay e 



^6,6 



(4.2) 
(4.3) 



Fx{XX)= [4,a 0] WM 



F{XX. 



Fy{XX)^ [0 h.h\ UM 



F{X,Y: 



F{X,Y) = QFi^,R), 

and M — QGQ^^. Ia,a and Ib^t are a x a and 6x6 identity matrices respectively. 

There is some freedom in the choice of 11, as long as the first a entries of the 
diagonal of A have zero real components. We define 11 as 



n 



n. 







where n+ £ M™ •'" is defined such that 

- |5i(n+diag(A+)n;i),J > - |5R(n+diag(A+)n;i), 

if ii < J2- 



Equations ^Jj and ([2^ show that F(0, 0) = 0, which gives F (0) = 0. Fur- 
thermore DF{0, 0) = (where DF is the Jacobian of F), so DF (0) = 0. Therefore 
we can apply centre manifold theory as described in |5] to predict the behaviour of 



(3.1 1 by projecting a system defined on C^™ '^™ down onto a much lower dimensional 



space of size C°'°. 

4.2. Centre manifold theory. Consider the system 

^X^AxX + Fx{X,H{X,e)) 
at 



(4.4) 



where H : 



■^a+l 



satisfies 



DH{X, e) [AxX + Fx{X, H{X))] = AyH{X, e) + Fy (X, H{X, e)), 

i/(0,0)=0, 
DH{0,0) =0. 
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Ax is a diagonal matrix whose entries have zero real parts and Ay is a diagonal 
matrix whose entries have non-zero real components. We expect the value of a to be 
relatively small. 

In the analysis of the suspended system shown below we project down onto the 
subspace spanned by the eigenvectors corresponding to the eigenvalues stored in Ax ■ 
This projection is carried out by the Q = I1P~^ matrix. By considering the block 
structure of the P^^ matrix we see that X^ is some linear combination of <f>fe and Rk. 
In other words, analysis carried out on mode fc in the reduced system will, directly, 
tell us about the behaviour of mode k in the full system (By the full system we mean 
the system defined in (2.9) and ( |2.10[ )). Or put more crudely, there is no mixing of 
the modes as we move between the two different spaces. 



By Theorem 2, Section 2.4 of "S], the solutions of (4.2 ) and (4.3 1 are stable (asymp- 
totically stable) (unstable) if the solutions of system (4.4) are stable (asymptotically 
stable) (unstable). 

Using the above conditions to find H is, in general, very difficult, instead we find 
an approximation ^ : C""*"^ — >■ C*". Define 



(7V*)(X, e) = D'i'iX, e) [AxX + Fx{X, *(X, e))] 
-(Ay*(X,e) + Fy(X,*(X,e))). 



(4.5) 



Suppose *(0, 0) = 0, i:)*(0, 0) = and that 



|(iV*)(X,e)||, = 



as II \X'^ e] ||q+i — > where g > 1. From Theorem 3, Section 2.5 of [5], 

\\H{X,e)-^{X,e)\\, = 



||.||a+i and ||.||f, are some norms defined on C""*"^ and C** respectively. 

For now, lets assume ^(X, e) = e'^LX + '^N{X) where ^l S C*^" is a matrix that 
will be determined below, and ^jv : C^ ^- C° is a quadratic of the form '^ i^{X) = 
[0 hb] Tn{X) where (TAr(X))fc - E,-, E,-, Cf,-,j,)^.-i^.V Rewrite M - QGQ-^ 
as 



M = 



M21 M22 



where Afn € C^" . 

To ensure a good approximation to the centre manifold H we want to make 



(7V*)(X,e) in Equation (4.5) 'small'. After substituting *(X,e) into (iV*)(X,e), 



we firstly collect all of the terms that are linear in X. These terms are 

e^L^xX + e^^lA/iiX + e3*LA/i2*LX - (eAy^iX + eKhiX + e^M22'^LX) . 

As Ax, Ay, Mil and M22 are all diagonal matrices it is straight forward to find ^^ 
so that it solves the equation 



OHlAxX + e^^iMiiX - (eAy*L^ + eK'hiX + e^M22'^LX) 



0. 



This would leave us with terms that are 0{e^X). However, to help simplify the 
analysis of the reduced system we instead choose ^P^ so that it solves e'i'j^AxX — 
{sAy^lX + eM2iX) = 0. In particular, we set 



(*l). 



(M21), 



''' {kx)j,-{KY)u 

The remaining terms are 0{€^X). 

Referring back to {N'^){X , e), the terms that are second order in X are 

e2*iAfi2*w(X) + e*L [h^a 0] F(X, e^^X) 

+ D^N{X)kxX + eL'*Ar(X)AfiiX 

+ e^D'^N{X)Mi2'^LX 

- (Ay*Ar(X) + eAf22*w(X) + [0 h.b] F{X, s^lX)) 

Setting Tn{X) so that 



DTNiX){Ax+eM,,)X- 





(Ay + eM22) 

leaves remainder terms that are 0{e'^XX'^). Now 



Tx{X) = F{X,e'^LX), (4.6) 



(DTxiX)),^^ = Y: [ef,,,,) 



khJ) 



Xj2- 



32 



Consider a mode k such that (Ay)^ is defined. The entry of the left hand side of 
Equation (4.6) corresponding to mode k can be written as 



zZ zZ ^{3u32Ajij2)^n^. 



32 J 



(4.7) 



^^^^^ 4i J2) ^ (^^ + ^^'^ii)(ji,ji) + (^x + ^^'^ii)(j2,J2) - (^^^ + f^'^22)(fe,fe). Using 



= Q- 



a.a 



with Equations (2.7) and ( |2.8[ ) we see that 

{F{X,e^,X))^^Y.J2fL32)XnXj 



(4.8) 



where the constants fh ■ s are known (can be calculated). Comparing Equations 
{ il\ and ( [Xs] ) suggests 

ffc _ /-fc I k 

'^{31,32) J (31,32)1 {31, 32)- 

So ^(X, e) = €'^ lX + \E'7v(-X') gives an third order approximation to H. 
Substituting ^ back into (4.4), the reduced system becomes 

d 



dt 



X - A^X + e (AfnX + Mi2*(X, e)) + [/,,, O] F (X, *(X, e)) 



= [AxX + eAfnX + e^Mi2-^LX] 

+eA-h2^x{X,e)+[Ua O] F (X, ^(X, e)) . (4.9) 

We can study the behaviour of the above a x a system to infer the behaviour of the 
original system in (4.2) and (4.3), and consequently (2.9) and (2.10). 
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Fig. 5.1. Maximum real part of the non-zero eigenvalues of L as a is increased with different 
values of /3 



5. Application of centre manifold analysis - periodic boundaries. To 

help better understand the behaviour of the Hasegawa-Wakatani equations, the non- 
linear analysis developed in Section |4.2| is firstly applied to example problems with 
periodic boundary conditions. The zero boundary case is studied in Section [6) the 
results differ greatly from those presented in this section. 

In Section |5.1| the modes are analysed individually. That is, in the centre man- 
ifold analysis we only concentrate on those modes whose eigenvalues have zero real 
components for a given values of a. In Section [5.2| we present some experiments that 
are designed to take the interaction between the modes into account. 

We will ignore the fc — (0, 0) mode in the following discussion as it does not 
influence the behaviour of either the full or reduced systems. To simplify the discussion 
we set Pp = j3^ = (3. 

5.1. Individual mode study. In the first example we set n = 16, k = 1.5, p = 
1 and try varying a. The G matrix is as defined previously in Section [4. 1[ 

A small script was written in Scilaqjto evaluate the eigenvalues. The plots in 
Figure |5.1| show the maximum value of the real part of the non-zero eigenvalues of L 
for /3 = 0.1, 0.05, 0.01, 0.005 and 0.001, and varying a. In the centre manifold theory 
we are interested in the choices of a that give eigenvalues with zero real components. 
It may not be so obvious in the plots, but the graphs cross the x-axis twice, once very 
close to the y-axis. That is, for each /3 there is an q;„ where the 5ft (A^ ) < for all fc 
if a > a„ and there is an a/ where the 5ft (A^ ) < for all fc if a < a;. In other words, 
the solutions converges to zero if a > a„ and a > a; . Consequently, we are interested 
in the region ai < a < a„. We concentrate on the results for larger values of a and 
now set /3 = 0.001. 

In this case L is a matrix of size 2046 x 2046. If we were unable to exploit the 
sparse structure of the matrix and had to calculate the eigenvalues and eigenvectors 
numerically the computational cost would be exorbitant. 



^ http : //www . scllab . org/ 
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The numerical simulation of the original full system in (2.9 1 and (2.10) were 
carried out using the Fortran 90 code on a grid of size 64 x 64. 

For the initial condition $", the coefficients of <i>^ were chosen randomly if < 
\kx\ < 5 and 1 < \ky\ < 5, otherwise <f>^ is zero. Furthermore |<i>^| < 0.01, i?" = ^° 

and $^,^,,^) = $?-fe,._fe,) (for aU k). 

Starting with a large value of a and slowly reducing it, we found the first eigenval- 
ues to have zero real components are At, _^^s when a = 281.2475. (So, a„ =281.2475.) 
According to the centre manifold analysis the reduced system is the following linear 
system (we are ignoring the k = (0, 0) mode) 

X(o,_i) = (7.500 X lO^^i + e(-3.556 x 10"*^ + 1.896 x IQ-^i)) X(o^-i), 



dt 



e^ (6.321 X 10"^ - 5.899 x 10""i) X(o,_i), (5.1) 






The first point to note is that the nonlinear terms in Equation (4.9) are zero for 
this particular example. To see why, firstly observe that the vector [$x Rx] = 
Q^^ \X 'i'{X,e)] will be zero everywhere except for those entries corresponding 
to the k = (0, ±1) modes. Therefore the inverse Fourier transform of $x and Rx are 
one dimensional, real, functions of y. The Poisson bracket applied to such functions 
is zero. 



If a ^ 281.2475 (e ^ 0), the real part of the eigenvalues of Equation (5.1) are 
negative, indicating that the solution should converge towards zero. The numerical 
solution of both the reduced and full systems did converge towards zero. 



Given the small values of the real components of the eigenvalues of (5.1), the 



solutions of the full system should be stable when a ^ 281.2475. Figure 5.2 shows the 
real part of $(o,i) and R(o,i) when a = 280. Note however, as a is decreased further 
the real components of the eigenvalues will increase, and possibly become positive, so 
for certain values of a we may see exponential growth in the k = (0, ±1) modes. 

Continuing to reduce a we notice 3^? ( ^^iti.]^ _|_i-) ) = at a = 83.326665. The 
reduced system is 

X(_i _i) = (5.000 X IQ-^i + e(-2.400 x 10"^ + 3.840 x lO^'^i)) ^(-i,-i) 

+6^ (1.919 X lO"'^ - 5.375 x IQ-^i) X(_i,_i) 

+e (Mi2*A,(X, e))(_i,_i) + F{X, vI/(X, e))(_i,.i), (5.2) 



dt 



d^ 

dt"'-'''' Jt 



— X(_i_i) = — X(_i,_i), 






To solve Equation (5.2) we used the Scilab ordinary differential equation solver. 

For e > the eigenvalues of the linear part of the reduced system have nega- 
tive real components, so the solution should converge towards zero. The numerical 
solutions of the reduced system converged to zero for such values of e. 
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Fig. 5.2. Real part o/*(o,i) and K(o,i) M a = 280 and 0<t< 5000. 
reduced equation at epsilon = -0.5 
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Fig. 5.3. Real part of Xr^i ^i-. when e = —0.5 



The real part of ^(i.i) is shown in Figure 5.3 for e = —0.5 The resuhs imply 
the k — (±1,±1) modes are stable when a « 83.326665, which was the case in the 
numerical simulations of the full system. 

The fc = (±1,±1) mode may be stable, but according to Equation (5.1) the 
k = (0, ±1) mode will grow when a ~ 83.326665, as shown in Figure 5.4 



Reducing the value of a further the next eigenvalue to have real zero components 
at a = 71.984. As with the k — (0,±1) case, the nonlinear term drops out 
and the reduced system is 



^^ •^(0,±2) 






(5.999 X lO^^i + e(-5.554 x lO"'"^ + 1.482 x IQ-^i)) X(o,_2) 
+6^ (6.165 X 10"^ - 2.878 x lO^^z) X(o,_2), 



(5.3) 



dt 



:X. 



(0,2) 



dt 



(0,-2)- 



For negative values of e (a < 71.984) the eigenvalues of (5.4) have positive real 



components so we expect to see a growth in the fe = (0, ±2) modes. 



The next eigenvalue to obtain zero real components is A 



+ 

(±i,±2; 



at a = 41.64583. 



We used the Scilab's ordinary differential equation solver to find the solutions of the 
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Fig. 5.5. Real part o/^(±i,±2) when e = —0.5. 



reduced system when e < 0. The results in Figure [575] indicate this is a stable mode 
when a ~ 41.64583. Once again, exponential growth in the k ~ (0,±1) mode is 
evident when a = 41. 

It is important to carry out the numerical simulations on appropriate sized grids. 
For example, if we rerun the previous experiment on a grid of size 32 x 32, instead of 
64 X 64, the k = (0, ±1) mode appears to be stable, as shown in Figure 5.6 This is 



incorrect, the solution should grow exponentially. If we did not have the analytical 
results and just looked at the results for the fc — (0,±1) mode shown in Figure 5.6 



it would not be obvious that something is wrong. But, the apparent high frequency 



(noisy) solution of $(±i,±i) and i?(±i,±i) in Figure 5.7 is an indication of the presence 
of numerical errors. 

The above results show how important it is to use a computational grid of the 
correct size. If the grid is too small the qualitative type of information extracted from 
the solutions can be very misleading. 

Reducing a even further, the next eigenvalues to have zero real components are 
'^(o,±3) ^h'^n " = 20.20947, (±1, ±3) when a = 15.168627, (±2, ±2) at a = 12.310091, 
(±2, ±1) at a = 10.39582 etc. In other words, as a is reduced more modes influence 
the long term behaviour of the system. The size of the computational grid must 
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Fig. 5.7. Real part of ^i\^i), and Rn^\) for a = 41, < t < 5000 and numerical grid of size 
32 X 32, 



be big enough to ensure these modes are accurately represented in the numerical 
computations. 



5.2. Solution of a model problem. The experiments in Section 5.1 focussed 
on one particular mode at a time. Here we concentrate on a choice of parameters 
where a large number of modes are expected to contribute to the long term behaviour 
of the system. In particular, we set p = 2, a = 1 and k = 1. The experiments were 
run using different values of /3 (J3 — 10"^, (3 — 10^'^ and (3 = 10^^). We also tried 
computational grids of size 64 x 64 to 256 x 256 to ensure that results are consistent. 
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Fig. 5.8. Real part o/<I>(q ^j, and ^(o,i) /c" 7 = 10 ^ and < i < 350. 



We know, analytically, that if the initial conditions $* and R* are chosen so that 
$^ and i?^ are zero everywhere except the k = (0,±1) mode, the solution will grow 
exponentially. This is also evident in the numerical simulations. Is this behaviour 
restricted to very specific choices of initial conditions? Or is it an indication of what 
we should expect in general? 

Lets consider another set of initial conditions 4>° and RP where the coefficients 
of ^^ are chosen randomly if < \kx\ < 5 and 2 < \ky\ < 5; and 'fm ±i) ~ ^•-'^■ 
Otherwise $^ is zero. Furthermore, the values are scaled if A; ^ (0,±1) so that 
I'I'fel < 7 ■ For small values of 7 the initial condition is close to $* and R* . Finally 

iJ" = *Oand$«,^,,^)=$?_,^,_,^). 

We found that for 7 := 10 ® and 7 = 10 ^ the solution showed exponential growth 
for all values of /?. In other words, any initial conditions close to $* and R* appear 
to be unstable. 

When 7 = lO^'^ the solution once again showed exponential growth for f3 = 10^^, 
but it behaved differently when f3 < 10"'^. As shown in Figure 5.8 the k = (0,±1) 
modes initially appeared to be stable, until about t = 300 when they once again grow 
exponentially. The growth appears to be primarily in the k — (0, ±1) modes as no 
growth is was evident in the k = (±1, ±1) modes. 

We also tried other initial conditions and obtained similar results to those men- 
tioned above. We thus conclude that for periodic boundary conditions, the nodes 
sitting on the (0, ±fcy) boundary will always grow for certain choices of a. 

6. Application of centre manifold analysis - zero boundary conditions. 

In Section [5] we considered the case of periodic boundary conditions, in this section 
we assume that the y boundaries are zero and the x boundaries are periodic. Zero 
boundary conditions can be enforced by the use of sin transformations, which may 
in-turn be implemented through Fourier transformations, see [2 . To implement a sin 
transformation of size N we need to use a Fourier grid of size 2N. 

When zero boundary conditions are applied along the y-axis the k = (0, ky) modes 
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Fig. 6.1. Real part o/#(i ij and -R(i,i) for a = 70 and <t < 5000. 



will always be zero. We saw in Section [5] that these modes played a special role. We 
expect to see a marked changed in the behaviour of the system when zero boundary 
conditions are enforced. 

To find the initial conditions we firstly evaluate 



P 



4 m-1 

EE 

k=l a,b=0 



^kVk si'n{2TT ka/m) su\{2Tikh/m) 



ni—l 

E 

a=0 



cos{27rka/m) 



where 7fc and rjk are chosen randomly. We then set $ = sJ^^^p/||J^^^p||i where J- is 
the Fourier transformation, ||.||i is the LI norm and s is a scaling term. Furthermore 
i? = $ and $fi. I. -I = ^(-h -k )■ 

6.1. Individual mode study. In this section we use the same parameters as 
in Section |5.1[ but apply zero boundary conditions along the y-axis. The numerical 
calculations were carried out on a grid of size 128 x 128. The scaling parameter s is 
0.01. 

When a > 83.326665 the real part of all of the eigenvalues of L are negative, so 
the solution should converge to zero. The numerical simulations converged to zero. 

We know that 3? [^t±i +i) ) = at a = 83.326665. The solution of the reduced 
system shown in Figure 5.3 suggests the k — (±1,±1) modes are stable when a « 
83.326665. The modes were stable in the numerical simulations. 



However, we see from Equation (5.2) that as a is decreased the real part of the 



eigenvalues of the linear system increase, and as shown in Figure [6T] this growth may 
dominate the system. Such growth in the k = (±1,±1) modes was also confirmed by 
numerical solutions of the reduced system. 

Although the results for a = 70 suggests the k = (±1,±1) will always grow 
exponentially, the plot for a = 41 shown in Figure [6?2] contradicts such an assumption. 
We do not believe that the change in behaviour at about t = 3000 is a result of 
numerical error, but rather an example of the stable bifurcation described in Section 
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Fig. 6.2. Real part o/f&ji jj and -R(i,i) for a = 41 and <t < 5000. 
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Fig. 6.3. Real part of^,-j^m and Rnfi) for a = 41 o,nd <t < 5000. 



|6.2[ It was necessary to increase the size of computational domain to 512 x 512 to 
accurately capture the results as many of the higher order modes grew in the region 
where the system suddenly changed its behaviour. For relevance to the discussion in 



Section 6.2 we have also included a plot of the k ~ (±1, 0) modes in Figure 6.3 



6.2. Solution of a model problem. The model problem under consideration 
here is the same as the one in Section[5?2l but with the y-axis set to zero. Consequently, 



the (0, ±1) modes that showed exponential growth in Section 5.2 have been removed 



from the system. The initial conditions are the same as those mentioned in the start 
of Section [6l 
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Fig. 6.5. Real part o/$(i.i) and -R{i,i) for /3 = 0.001 and <t < 3000. 



The properties of the solution varied greatly with /3. Figure |6.4| shows the real 
part of $(1,1) and -R(i,i) for /3 = 0.01. According to the solution of the corresponding 
reduced system, the k = (±1, ±1) mode should be stable, which is in agreement with 



the numerical results shown in Figure 6.4 



Figure 6.5 shows the real part of $(i,i) and -R(i,i) for (3 = 0.001. This time the 
solution of the reduced system shows a growth in the k = (±1,±1) modes. We see 



from Figures 6.5 and 6.6 that the (±1, ±1) and (±1, 0) modes interact in such as way 



as to result in stable bifurcation. 
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Fig. 6.6. Real part o/'I>(i q) and R(ifi) for f3 = 0.001 and < t < 3000. 



7. Conclusion. Our initial motivation for using the centre manifold theory to 
study the nonlinear behaviour of Hasegawa-Wakatani equations was to verify the 
results of our code. The theory proved to be a useful tool in predicting and explaining 
the behaviour of the results. It also highlighted some unexpected properties of the 
Hasegawa-Wakatani equations. 

We showed that it is important to carry out the computations on the correct 
sized grids, otherwise the qualitative behaviour of the solution will be wrong. One 
advantage of working in the Fourier space is that it is easy to see when a solution 
is wrong. A noisy solution like the one shown in Figure |5.7| is a strong indication of 
numerical errors. The author is not aware of any study that has been carried out 
to determine what are appropriate sized grids for finite element or finite difference 
discretisations. 

We presented a number of different examples of unstable solutions that arose 
when using periodic boundary conditions. Some of our results directly contradict the 
stable solutions reported in other papers. The experiments presented here suggest 
that any numerical simulation showing stable solutions should be checked to verify 
they are not a consequence of numerical errors or a consequence of terminating the 
simulations too soon. 

We showed examples of stable bifurcation behaviour occurring when the boundary 
along the y-axis was set to zero. Our future work will focus on better understand- 
ing how the different parameters influence the bifurcation behaviour. We intend to 
carry out this analysis by once again using the centre manifold analysis, but this 
time projecting down onto larger subspaces containing more modes of interest, as 
well as continuing to verify our results by numerical solutions on the full system of 
equations. 
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